#pragma once
#include <cstdio>
#include <functional>
#include <utility>
#include <cmath>
#include <sys/cdefs.h>
#include "utils.hpp"
#include "treduce.hpp"
namespace seqfilt {
template <int I, int FO>
struct filter_out_t {
  typedef utils::iseq<int, I> out_seq;
};
template <int I>
struct filter_out_t<I, I> {
  typedef utils::iseq<int> out_seq;
};
template <int I, int F>
using filter_out_one = typename filter_out_t<I, F>::out_seq;

template <int... Is>
auto merge_seq(utils::iseq<int, Is...> seq0) {
  return seq0;
}
template <int... Is, int... Js, typename... Ts>
auto merge_seq(utils::iseq<int, Is...> seq0, utils::iseq<int, Js...> seq1, Ts... rest) {
  auto ijseq = utils::iseq<int, Is..., Js...>();
  return merge_seq(ijseq, rest...);
}
template <int F, int... Is>
auto filter_out(utils::iseq<int, Is...> seq) {
  return merge_seq(filter_out_one<Is, F>{}...);
}
}; // namespace filterseq

template <typename T, int N, int M>
struct matrix;
template <bool F2D, bool N1>
struct wrap_initializer {
  template <class F>
  static auto wrap(F f) {
    return f;
  }
};
template <bool N1>
struct wrap_initializer<false, N1> {
  template <class F>
  static auto wrap(F f) {
    return [&](int i, int j) { return N1 ? f(j) : f(i); };
  }
};
constexpr int const_ctz(int c) {
  return __builtin_popcount((c & -c) - 1);
}

template <typename T, int N>
using colvec = matrix<T, N, 1>;
template <typename T, int N>
using rowvec = matrix<T, 1, N>;

template<int N, typename T>
__always_inline T matdet(T *mat, utils::iseq<int> iseq, utils::iseq<int> jseq) {
  return 1;
}

template <int N, typename T, int I, int... Is, int... Js>
__always_inline T matdet(T *mat, utils::iseq<int, I, Is...> iseq, utils::iseq<int, Js...> jseq) {
  auto iseq_next = utils::iseq<int, Is...>{};
  auto kseq = jseq;
  return treduce::addsub(mat[I * N + Js] * matdet<N>(mat, iseq_next, seqfilt::filter_out<Js>(kseq))...);
}

template <int N, typename T>
__always_inline T matdet(T *mat) {
  auto nseq = utils::make_iseq<int, N>();
  return matdet<N>(mat, nseq, nseq);
}
template <typename... Ts>
static __always_inline void eval_all(Ts... args) {
}
template<int N, int IJ>
static constexpr int cosign() {
  constexpr int i = IJ / N, j = IJ % N;
  return (i + j & 1) ? -1: 1;
}
template <int N, typename T, int ...IJs>
__always_inline void matinv(T *ret, T *mat, utils::iseq<int, IJs...> ijseq) {
  auto iseq = utils::make_iseq<int, N>();
  T value1d[] = {cosign<N, IJs>()*matdet<N>(mat, seqfilt::filter_out<IJs%N>(iseq), seqfilt::filter_out<IJs/N>(iseq))...};
  T determ = matdet<N>(mat);
  if ((T)0.5 == 0.5) {
    T invdet = 1 / determ;
    eval_all(ret[IJs] = value1d[IJs] * invdet...);
  } else {
    eval_all(ret[IJs] = value1d[IJs] / determ...);
  }
  // return matrix<T, N, M>([&](int i, int j) {return value1d[j * N + i] / determ;});
}
template<int N, typename T>
__always_inline void matinv(T *ret, T *mat) {
  auto ijseq = utils::make_iseq<int, N*N>();
  return matinv<N>(ret, mat, ijseq);
}

template <int ...Js, typename T>
__always_inline T tdot_expand(T *v0, T *v1, utils::iseq<int, Js...> jseq) {
  return treduce::sum(v0[Js] * v1[Js]...);
}
template <int ...Js, typename T>
__always_inline T tdot_seq(T *v0, T *v1) {
  return treduce::sum(v0[Js] * v1[Js]...);
}
template <int M, typename T>
__always_inline T tdot(const T *v0, const T *v1) {
  return tdot_expand(v0, v1, utils::make_iseq<int, M>{});
}
#undef min
#undef max
#ifdef __sw_host__
#include <iostream>
template<int N, typename T>
void printmat(T *mat) {
  for (int i = 0; i < N; i ++) {
    for (int j = 0; j < N; j ++) {
      std::cout << mat[i * N + j] << " ";
    }
    std::cout << std::endl;
  }
}
#endif
#ifdef TEST
#include <cstdio>

int main(){
  double d[] = {6, 1, 1, 4, -2, 5, 2, 8, 7};
  printf("%d\n", matdet<3>(d));
  double dinv[9];
  matinv<3>(dinv, d);
  printmat<3>(dinv);
  double mult[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
  for (int i = 0; i < 3; i ++) {
    for (int j = 0; j < 3; j ++){
      for (int k = 0; k < 3; k ++) {
        mult[i * 3 + k] += dinv[i * 3 + j] * d[j * 3 + k];
      }
    }
  }
  printmat<3>(mult);
}
#endif
// template <int ...Is, int ...Js, int M, typename T>
// void tmatvec_expand(T *ret, T (*mat)[M], T *vec, utils::iseq<int Js...> jseq) {
//   ret[Is]
// }
// template <int N, int M, typename T>
// void tmatvec(T *ret, T (*mat)[M], T *vec) {
//   // auto iseq = utils::make_iseq<int, N>();
//   auto jseq = utils::make_iseq<int, M>();
// }
// int main(int argc, char **argv) {
//   auto f = [](int i, int j) { return rand()%16; };

//   matrix<double, 4, 4> a(f);
//   // auto g = [](int i, int j) { return i * 2; };
//   // auto h = [](int i) { return i * 2; };
//   // colvec<int, 4> b(h);
//   // auto c = a.matmul(b);
//   matrix<double, 3, 3> d;
//   // d(0, 0) = 6;
//   // d(0, 1) = 1;
//   // d(0, 2) = 1;
//   // d(1, 0) = 4;
//   // d(1, 1) = -2;
//   // d(1, 2) = 5;
//   // d(2, 0) = 2;
//   // d(2, 1) = 8;
//   // d(2, 2) = 7;
//   d(0, 0) = 3;
//   d(0, 1) = 0;
//   d(0, 2) = 2;
//   d(1, 0) = 2;
//   d(1, 1) = 0;
//   d(1, 2) = -2;
//   d(2, 0) = 0;
//   d(2, 1) = 1;
//   d(2, 2) = 1;
//   auto e = d.inv();
//   // printf("%d\n", treduce::addsub(1, 2, 3, 4));
//   printf("%d %d\n", d.det(), a.det());
//   // printf("%d\n", d.tdet(utils::make_iseq<int, 3>{}, utils::make_iseq<int, 3>{}));
// }
// int main(int argc, char **argv) {
//     auto f = [](int i, int j) {return i + j;};
//     int a[4][4];
//     make_matrix<4, 4>(a, f);
//     for (int i = 0; i < 4; i ++) {
//       printf("[");
//       for (int j = 0; j < 4; j ++) {
//         printf("%d, ", a[i][j]);
//       }
//       puts("],");
//     }
//     auto g = [](int i) {return i*2;};
//     int b[4];
//     make_array<4>(b, g);
//     for (int i = 0; i < 4; i ++) printf("%d\n", b[i]);
//     int c[4];
//     matvec<4, 4>(c, a, b);
//     for (int i = 0; i < 4; i ++) printf("%d\n", c[i]);
//     return 0;
// }
// template <int I, int J, class F, class T, int N, int M>
// int make_element(matrix<T, N, M> mat, F f) {
//   mat[I][J] = f(i, j);
//   return 0;
// }
// template <int I, int J, class F, class T int N, int M>
// int make_row
// template <int I, int J, int NJ, class F>
// int make_element(std::result_of_t<F(int, int)> ret[][NJ], F f) {
//     ret[I][J] = f(I, J);
//     return 0;
// }
// template <int I, class F>
// int make_element(std::result_of_t<F(int)> ret[], F f) {
//     ret[I] = f(I);
//     return 0;
// }
// template<typename ...Ts>
// void runall(Ts... args){

// }
// template <int I, int ...Js, class F>
// int make_row(std::result_of_t<F(int, int)> ret[][sizeof...(Js)], F f, utils::iseq<int, Js...> jseq) {
//     runall(make_element<I, Js>(ret, f)...);
//     return 0;
// }
// template <int ...Is, int ...Js, class F>
// void make_rows(std::result_of_t<F(int, int)> ret[][sizeof...(Js)], F f, utils::iseq<int, Is...> iseq, utils::iseq<int, Js...> jseq) {
//     runall(make_row<Is>(ret, f, jseq)...);
// }
// template<int NI, int NJ, class F>
// void make_matrix(std::result_of_t<F(int, int)> ret[][NJ], F f) {
//     auto iseq = utils::make_iseq<int, NI>{};
//     auto jseq = utils::make_iseq<int, NJ>{};
//     make_rows(ret, f, iseq, jseq);
// }
// template<int ...Is, class F>
// void make_elements(std::result_of_t<F(int)> ret[], F f, utils::iseq<int, Is...> iseq){
//   runall(make_element<Is>(ret, f)...);
// }
// template<int NI, class F>
// void make_array(std::result_of_t<F(int)> ret[], F f) {
//   auto iseq = utils::make_iseq<int, NI>{};
//   make_elements(ret, f, iseq);
// }

// template<typename T0, typename T1>
// struct add_result{
//   typedef decltype((T0)0 + (T1)0) type;
// };
// template<typename T0, typename T1>
// using add_type = typename add_result<T0, T1>::type;
// template<typename T0, typename ...Ts>
// struct sum_result{
//   typedef add_type<T0, typename sum_result<Ts...>::type> type;
// };
// template<typename T0>
// struct sum_result<T0> {
//   typedef T0 type;
// };
// template<typename ...Ts>
// using sum_type = typename sum_result<Ts...>::type;
// template<typename T>
// T sum(T first) {
//   return first;
// }
// template<typename T, typename... Ts>
// sum_type<T, Ts...> sum(T first, Ts... args) {
//   return first + sum(args...);
// }

// template<int I, int ...Js, typename T>
// int matvec_element(T ret[], T mat[][sizeof...(Js)], T vec[], utils::iseq<int, Js...> jseq){
//   ret[I] = sum<int, int, int, int>((mat[I][Js] * vec[Js])...);
//   return 0;
// }
// template<int ...Is, int ...Js, typename T>
// void matvec_expand(T ret[], T mat[][sizeof...(Js)], T vec[], utils::iseq<int, Is...> iseq, utils::iseq<int, Js...> jseq) {
//   runall(matvec_element<Is>(ret, mat, vec, jseq)...);

// }
// template<int NI, int NJ, typename T>
// void matvec(T ret[], T mat[][NJ], T vec[]) {
//   auto iseq = utils::make_iseq<int, NI>{};
//   auto jseq = utils::make_iseq<int, NJ>{};
//   matvec_expand(ret, mat, vec, iseq, jseq);
// }
